home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / fft.lha / fft / fft.h < prev    next >
C/C++ Source or Header  |  1993-08-08  |  5KB  |  156 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *                  Fast Fourier Transform
  6.  *
  7.  * The present package is intended to be used for computing the sums
  8.  * of the form
  9.  *    (1) xf[k] = SUM{ x[j] *  exp(-2*PI*I/N j*k), j=0..N-1 }, k=0..N-1
  10.  *    (2) x[j]  = 1/N SUM{ xf[k] * exp(+2*PI*I/N j*k), k=0..N-1 }, j=0..N-1
  11.  *    (3) xf[k] = SUM{ x[j] *  exp(-2*PI*I/N j*k), j=0..N/2-1 }, k=0..N/2-1
  12.  *
  13.  * where N is the exact power of two.
  14.  *
  15.  * Formula (1) defines the classical Discrete Fourier transform, with x[j]
  16.  * being a real or complex sequence. The result xf[k] is always a complex
  17.  * sequence, though the package lets user to get only real and/or imaginaire
  18.  * parts of the result (sine/cosine transform of x if x is real), or the
  19.  * absolute value of xf (the power spectrum).
  20.  * Formula (2) is the inverse formula for the DFT.
  21.  * Formula (3) is nothing but the trapezoid rule approximation to the
  22.  * Fourier integral; the trapezoid rule is the most stable one with respect
  23.  * to the noise in the input data. Again, x[j] maybe either real or
  24.  * complex sequence. The mesh size other than 1 can be specified as well.
  25.  *
  26.  ************************************************************************
  27.  */
  28.  
  29. #pragma once
  30. #pragma interface
  31.  
  32. #include "LinAlg.h"
  33. #include <math.h>
  34. #include <Complex.h>
  35.  
  36. class FFT
  37. {
  38.   const int N;                // No of points the FFT packet
  39.                     // has been initialized for
  40.   int logN;                // log2(N)
  41.   const double dr;            // Mesh size in the r-space
  42.  
  43.   Complex * A;                // [0:N-1] work array
  44.                     // the transform is placed to
  45.  
  46.   Complex * A_end;                 // Ptr to the memory location next to
  47.                     // the last A element
  48.  
  49.   short * index_conversion_table;    // index_conversion_table[i]
  50.                     // is the bit-inverted i, i=0..N
  51.  
  52.   Complex * W;                // FFT weight factors
  53.                     // exp( -I j 2pi/N ), j=0..N-1
  54.  
  55.             // Private package procedures
  56.   void fill_in_index_conversion_table(void);
  57.   void fill_in_W();
  58.  
  59.   void complete_transform(void);
  60.  
  61. public:
  62.  
  63.   FFT(const int n, const double dr=1);    // Constructor, n being the no. of
  64.                     // points to transform, dr being the
  65.                     // grid mesh in the r-space
  66.   ~FFT(void);
  67.  
  68.             // Fundamental procedures,
  69.             // Input the data and perform the transform
  70.   void input(                // Preprocess the real input sequence
  71.          const Vector& x);        // Real [0:N-1] vector
  72.  
  73.   void input(                // Preprocess the complex input seq
  74.          const Vector& x_re,    // [0:N-1] vector - Re part of input
  75.          const Vector& x_im);    // [0:N-1] vector - Im part of input
  76.  
  77.             // Preprocess the input with zero padding
  78.   void input_pad0(            // Preprocess the real input sequence
  79.          const Vector& x);        // Real [0:N/2-1] vector
  80.  
  81.   void input_pad0(            // Preprocess the complex input seq
  82.          const Vector& x_re,    // [0:N/2-1] vector - Re part of input
  83.          const Vector& x_im);    // [0:N/2-1] vector - Im part of input
  84.  
  85.  
  86.             // Output results in the form the user wants them
  87.   void real(                // Give only the Re part of the result
  88.         Vector& xf_re);        // [0:N-1] vector
  89.  
  90.   void imag(                // Give only the Im part of the result
  91.         Vector& xf_im);        // [0:N-1] vector
  92.   
  93.   void abs(                // Give only the abs value
  94.        Vector& xf_abs);        // [0:N-1] vector (power spectrum)
  95.  
  96.                 // Return only the half of the result
  97.                 // (if the second half is unnecessary due
  98.                 // to the symmetry)
  99.   void real_half(            // Give only the Re part of the result
  100.          Vector& xf_re);    // [0:N/2-1] vector
  101.  
  102.   void imag_half(            // Give only the Im part of the result
  103.          Vector& xf_im);    // [0:N/2-1] vector
  104.   
  105.   void abs_half(            // Give only the abs value
  106.          Vector& xf_abs);    // [0:N/2-1] vector
  107.  
  108.  
  109.             // Perform sin/cos transforms of f: R+ -> R
  110.             // Source and destination arguments of the functions
  111.             // below may point to the same vector (in that case,
  112.             // transform is computed inplace)
  113.  
  114.                          // Sine-transform of the function f(x)
  115.                 // Integrate[ f(x) sin(kx) dx], x=0..Infinity
  116.   void sin_transform(        // j=0..n-1, n=N/2
  117.     Vector& dest,            // F(k) tabulated at kj = j*dk
  118.     const Vector& src        // f(x) tabulated at xj = j*dr
  119.             );
  120.  
  121.                          // Cosine-transform of the function f(x)
  122.                 // Integrate[ f(x) cos(kx) dx], x=0..Infinity
  123.   void cos_transform(        // j=0..n-1, n=N/2
  124.     Vector& dest,            // F(k) tabulated at kj = j*dk
  125.     const Vector& src        // f(x) tabulated at xj = j*dr
  126.             );
  127.  
  128.                          // Inverse sine-transform of the function F(k)
  129.                 // 2/pi Integrate[ F(k) sin(kx) dk], k=0..Inf
  130.   void sin_inv_transform(    // j=0..n-1, n=N/2
  131.     Vector& dest,            // f(x) tabulated at xj = j*dr
  132.     const Vector& src        // F(k) tabulated at kj = j*dk
  133.             );
  134.  
  135.                          // Inverse cosine-transform of function F(k)
  136.                 // 2/pi Integrate[ F(k) cos(kx) dk], k=0..Inf
  137.   void cos_inv_transform(    // j=0..n-1, n=N/2
  138.     Vector& dest,            // f(x) tabulated at xj = j*dr
  139.     const Vector& src        // F(k) tabulated at kj = j*dk
  140.             );
  141.  
  142.  
  143.                          // Inquires
  144.   int q_N(void) const            { return N; }
  145.   int q_logN(void) const        { return logN; }
  146.   double q_dr(void) const        { return dr;  }
  147.   double q_dk(void) const        { return 2*PI/N/dr; }
  148.   double q_r_cutoff(void) const        { return N/2 * dr; }
  149.   double q_k_cutoff(void) const        { return PI/dr; }
  150.  
  151. };
  152.  
  153.  
  154.  
  155.  
  156.